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ACDOtJOOGMjT-  ABSTRACT 

This  thesis  proposes  a  method  of  identifying  the  dynamics  of  ship 
angular  motion  at  sea  as  the  basis  for  ship  motion  estimation.    Various 
sources  of  information  are  discussed.    A  digital  model  to  simulate 
ship's  motion  is  developed.    Identification  algorithms  are  presented 
with  the  results  of  their  stochastic  digital  simulation.     Sensitivity  of 
identification  error  to  sample  rate  was  investigated.    A  plan  for  ship- 
board implementation  utilizing  a  digital  computer  is  presented. 
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INTRODUCTION 

The  sophistication  of  modern  shipboard  systems  is  imposing 
accuracy  requirements  in  measuring  a  ship's  attitude.    A  10  mil  error 
in  yaw  correction  to  a  gun  platform  is  a  200  yard  error  at  10  nautical 
miles.    System  time  delays  and  additive  measurement  noise  are  key 
sources  of  system  error.    Automatic  carrier  landing  systems  require 
that  the  predicted  touch  down  point  on  the  deck  be  known.     Present 
system  effectiveness  degenerates  rapidly  as  sea  states  increase. 
When  the  wave-off  decision  point  is  reached  by  the  aircraft,  the  auto- 
mated landing  system  must  know  the  aircraft  position  relative  to  the 
desired  impact  point  at  the  time  of  touchdown.    If  pitch  and  yaw  can 
accurately  be  predicted,  a  true  all  weather  capability  may  be  antici- 
pated.   For  these  and  many  other  applications  the  accurate  knowledge 
of  present  position  and  the  prediction  of  future  positions  are  key 
elements  of  the  system  environment. 

An  accurate  prediction  of  ship  motion  must  necessarily  involve 
information  from  ±he  forcing  function.     Proper  interpretation  should  lead 
to  an  identification  of  the  current  dynamic  structure  of  the  ocean  waves, 
This  could  yield  a  more  accurate  input  to  the  prediction  of  sea  states 
and  ocean  weather  propagation. 

The  French  have  done  considerable  work  in  the  field  of  spectro- 
angular  ocean  wave  analysis  [1,2],     However  their  work  is  success- 
ful only  in  predicting  mean  and  extreme  wave  amplitudes  and  periods. 
The  primary  available  inputs  are  local  surface  wind  conditions  dis- 
tributed over  the  ocean  area .    Work  in  this  area  is  not  yet  sufficiently 
advanced  to  be  applied  to  a  real  time  estimation  problem, 

The  investigation  of  the  statistical  properties  of  ocean  waves  has 
been  a  well  traveled  route.    R.  L.  Wiegel  [  3  ]   from  the  University 
of  California  has  done  extensive  work  using  both  surface  and  bottom 
moored  pressure  sensing  units.    Data  reduction  again  shows  normal 
distributions  of  wave  period  and  amplitude  in  open  waters. 


Ship  motion  estimation  may  be  approached  as  a  stochastic  con- 
trol problem.    Thus  the  ship  is  the  plant  and  the  ocean  wave  is  the 
forcing  function  with  normal  Gaussian  distribution.     Pitch,  roll, 
and  yaw  will  be  treated  as  separate,  uncoupled  subsystems  dis- 
cretely sampled  to  facilitate  digital  simulation  and  data  processing. 
Only  the  roll  mode  will  be  referred  to  in  the  remainder  of  this  paper 
on  the  assumption  that  the  other  modes  may  be  estimated  by  the  same 
techniques . 


PLANT  SIMULATION 

To  evaluate  a  technique  for  ship  motion  estimation,  it  is  first 

necessary  to  have  a  source  of  data.    A  fourth  order  plant  composed  of 

two  cascaded  identical  second  order  systems  was  selected  for  motion 

simulation.     This  plant  has  known  natural  frequency  (  u)    )  and  damping 

n 

coefficient  (  C  ) , 


G(S)  T~2 " 2 

\s      +  2  f  w      s  +  (jQ 
\  n  n 

The  plant  may  be  represented  by 

X    =   AX  +  BW     where 


A  = 


1 


2  2 

s     +2£u)     s+u> 
n  n 


0 

1 

0 

0 

0 

0 

0 

1 

0 

and  B  = 

0 

0 

0 

0 

1 

0 

2 

-4  r 

3 
n 

-4 

2 
n 

(r2+ 

.5) 

-4  rw 

n 

1 

(1) 


(2) 


Selection  of  C  and  o>    determines  the  dominant  characteristics.    The 

n 


discrete  recursive  state  equation  in  matrix  form  is 
X(k  +  1)  =  0  X(k)  +  TW(k)  . 


(3) 


W(k)  is  the  discrete,  white  Gaussian  excitation  with  zero  mean.    The 
desired  system  output  is  the  discrete  time  sequence  of  roll  angle 
corresponding  to  X    where 

X^k)    =    HX(k)   ,     H  =  (1  0  0  ...  ). 
Additive  measurement  noise,  v(k),  may  be  included  at  the  output  to 
complete  the  simulation  model,  fig.   1.    V(k)  is  assumed  to  be  white 
Gaussian  noise,  uncorrelated  with  W(k) .    The  process  is  described  by 
z(k)  =  HX(k)  +  v(k).  (4) 

It  will  be  convenient  to  transform  0  into  the  canonic  companion 
matrix  form,  0  *,  using  Lee's  observability  criteria  [4  ]  ,  giving 


H 
H  0 


D  = 


H  0 


n-1 


*  = 


D  0  D 


-1 


0  : 
__j 

al    a2 


(5) 


and  the  companion  distribution  matrix 

r*  =  Dr, 


(6) 


From  the  work  of  Lee  and  Ho  [  5  ]   this  defines  a  new  recursion  equation 
with  the  same  output  values  as  equation  (4)  as  seen  below 


Y(k+1)  =  0  *  Y(k)  +  T  *  W(k) 
z(k)       =  H  Y(k)  +  v(k) 


(7) 


where  W(k),  H,  v(k),  and  z(k)  are  the  same  but  the  Y  elements  are 

different. 

Representing  the  simulation  model  in  this  form  reduces  both  the 

number  of  computations  and  the  time  required  to  generate  a  sample  set 

of  data  on  a  digital  computer.     Since  the  data  set  is  unchanged  by  the 

foregoing  manipulations,  the  dominant  characteristics  are  preserved. 

Values  of  T=  0.7  and  o>     =  0.2  77  were  used  as  reasonable  estimates 

n 

for  the  roll  dynamics  of  a  moderate  size  ship.     The  excitation,  W(k) , 
was  a  random  number  generating  subroutine  with  zero  mean  and  stan- 
dard deviation  of  one.     Multiplying  the  random  numbers  by  a  con- 
stant permitted  simulation  of  any  desired  sea  state.    An  example  of 
roll  angles  generated  using  a  constant  of  twenty  and  a  sample  period 
of  0 . 1  second  is  shown  in  fig.  2. 
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Fig.   1.     Block  Diagram  of  the  Simulation  Model. 


XrSCALE  -  1.88E+89  UMJT3/INCH. 

^ -scale  -  i.8aE+aQ  units/inch. 


Fig.  2.    Roll  Angles  vs  Time  (10  sec/in)  Generated  by  the  Simulation 
Model  Excited  with  Discrete  Gaussian  Noise,  a=  20,  at  the 
Sampling  Interval,   0.1  sec.     Dominant  Period  =  10  sec, 
r  =  0.7. 
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IDENTIFICATION 

The  problem  of  identification  arises  when  the  set  of  differential 
equations  describing  the  system  dynamics  are  not  known,  uncertain, 
or  are  subject  to  periodic  or  continuous  change.    On  a  large  time 
scale  the  ship  dynamics  may  be  expected  to  change  as  fuel,  stores, 
ammunition,  etc.  vary  in  quantity  and  location.    Short  time  changes 
may  be  anticipated  due  to  changes  in  ships  heading  and  speed. 

In  state  vector  notation  the  recursion  equation  for  ship  roll 
motion  (uncoupled)  will  have  the  form; 

X(k+1)  =  *X(k)  +  rw(k)  (8) 

where  W(k)  is  the  discrete  Gaussian  excitation.    $  is  the  fixed 
state  transition  matrix  of  rank  N,  the  order  of  the  system.    The  term 
fixed  is  used  here  on  the  assumption  that  the  changes  in  dynamics 
mentioned  above  occur  at  such  a  sufficiently  slow  rate  that  $  may  be 
adequately  represented  as  time  invariant  over  a  finite  time  interval. 
Thus  $  must  be  periodically  identified  to  follow  ship  motion  within 
accuracy  specifications. 
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IDENTIFICATION  OF  THE  FREE  DYNAMIC  SYSTEM 
A  free  dynamic  system  may  be  represented  by  the  equation 

Y(k+1)  =  0  *  Y(k)  .  (9) 

Lee  developed  an  identification  scheme  for  a  plant  having  no  numerator 
dynamics  [4  ]  .    Lee's  development  gave  the  following  equations; 


J  *  =  s   (s    y1 

zn     zn- 1 


(10) 


where 


2n-l 


z        z 

1        2 


z        z 

2         3 


n 


Zn     Zn+1       '  *    *  Z2n-1 


"Z2      Z3  •    • 

•  Vl" 

Z3      Z4 

• 

zn 

O                                   0 

• 

"                  .1                   •       •        •       • 

n+1 

'     "2n_ 

and  z(k)  =  HY(k) ,  a  scalar. 

Using  the  same  data,  the  order  of  the  system  may  also  be  iden- 
tified.    Choose  some  M  greater  than  N  and  build  an  M  x  L  matrix,  A. 
L  =  1,2,3 N,N+1, 


Zl  Z2        Zi 

Z2  Z3        Zl+l 

•  «  • 

•  •  • 

z  z  z 

m        m+1  m+Z-1 


A  = 


The  product  A  A  will  be  positive  definite  for  L  less  than  or  equal  to  N 
and  singular  for  L  greater  than  N.     Hence  the  order  of  the  system  is 
the  maximum  nonsingular  value  of  L. 

A  problem  arises  from  the  possible  occurrence  of  numerator  dy- 
namics which  leads  to  residues  (errors)  in  the  identification.    In  z- 
transform  theory  the  N      row  of  0  *  corresponds  to  the  denominator 
coefficients  of  the  system  z-transform.    The  numerator  remains 
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unidentified.    Hence,  if  there  is  an  unknown  non-scalar  numerator  in 
the  real  plant  corresponding  to  a  nonzero  initial  condition,  errors  will 
result. 

Identification  using  this  method  was  investigated  with  excellent 
results  by  setting  an  initial  condition  for  roll  angle  and  then  releasing 
the  plant.    If  excitation  is  applied  while  measurements  are  being  taken, 
the  plant  is  no  longer  a  free  dynamic  system.     Identification  attempted 
with  random  excitation  applied  was  inaccurate  and  unreliable. 
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IDENTIFICATION  FOR  THE  STOCHASTIC  CASE 
Using  the  transformed  recursion  equation  in  the  new  state  space, 
Lt.  Ralph  Hudson's  unpublished  doctoral  notes  show  the  development 
of  a  straight  forward  identification  scheme  as  follows: 


Y(k+1)  =  0  *  Y(k)  +  w(k) 


(11) 


where 


Y(k)  = 


'k-n+1 


'k-1 


,  <*j  (k)  =  r  *w(k) , 


and  the  scalar,  W(k) ,  represents  the  white  Gaussian  excitation,    If 


equation  (11)  is  post  multiplied  by  Y(k)    ; 

Y(k+l)Y(k)T  =  0  *Y(k)Y(k)T  +  03  (k)Y(k)T 
Taking  the  expectation  of  both  sides; 

E[Y(k+l)Y(k)T]  =  0  *E[Y(k)Y(k)T]   +  E[u)(k)Y(k)T] 


(12) 


(13) 


Notice  that  the  E[  Y(k+l)Y(k)    ]   is  just  the  autocorrelation  function, 

T 
R(  t),  for  tau  equal  to  one.     Similarly,  E[Y(k)Y(k)    ]   is  the  auto- 
correlation for  tau  equal  to  zero.    Therefore  0  *  =  R(1)R(0) 

If  the  order  of  the  identification  is  M  and  M  £N,  the  true  system 
order,  the  inverse  will  exist.  This  is  implied  by  the  fact  that  an  N 
order  linear  system  may  be  defined  by  N  linear  independent  differ- 
ential equations.  Hence  any  M  ^  N  of  these  equations  are  also  in- 
dependent. If  M  >N  independence  is  lost  and  R(0)  of  rank  M  >N  is 
singular.  Thus  if  N  is  unknown,  it  may  be  identified  by  testing  for 
the  largest  non-singular  rank  of  R(0). 

For  application  to  discrete  linear  (Kalman)  filtering  it  is  also 
useful  to  note  that  the  Q  required  for  determining  the  gain  matrix  may 
be  concurrently  identified  as  follows; 


Q  =T*  E[W(k)W(k)T]r*T  =  E  [  w  (k)  U3(k)T] 


(14) 
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Hence 


Q  =E  [{Y(k+1)  -  0*Y(k)}  {Y(k+1)  -  0*Y(k)}T] 


which  simplifies  to 


Q  =R(0)  -  0  *R(1)T.  (15) 


This  scheme  was  programmed  in  Fortran  63  for  simulation  on  the 
CDC  1604  digital  computer  using  data  from  the  simulation  model. 
Batch  processing  and  recursive  identification  were  investigated. 

1.  Batch  processing:    A  sample  set  of  roll  angles  was  generated 
and  stored.    The  stored  data  was  then  batch  processed  to  form  the  auto- 
correlation functions,  R(0)  and  R(l),  for  the  set.    A  matrix  inversion 
routine  was  used  to  invert  R(0).    0*  was  then  identified  as  the  matrix 
product,  R(1)R(0)~    . 

If  r  *  was  restricted  to  a  single  non  zero  element,  identifi- 
cation of  the  fourth  order  plant  was  good  using  50  samples  and  ex- 
cellent for  5  00  samples.    However,  if  T  *  had  all  non  zero  elements, 
residues  due  to  numerator  dynamics  caused  errors  as  great  as  110% 
of  the  true  values  of  the  elements.    Increasing  the  sample  size  from 
500  to  900  showed  no  improvement. 

2.  Recursive  identification:    Hudson's  identification  scheme  may 
be  manipulated  into  the  following  recursion  equations  by  using  a  matrix 
inversion  lemma. 

0*(k+l)  =  0*(k)  +  hf(k)TP(k)Y(k)  +  l  J  "1  fy(k+l)  -  0A*(k)Y(k)j  Y(k)TP(k) 
P(k+1)     =  P(k)  (  I  -    [  Y(k)TP(k)Y(k)  +  1  ]-1Y(k)Y(k)TP(k) 


where  P(k)  =    R(0)      . 


^Y(k)Tp(k)Y(k)  +  y ,~1, 


Let  j8  (k)  =   (Y(k)   P(k)Y(k)  +  1  )       ,  a  scalar. 

Then 

0*(k+l)  =  0A*(k)  +  |8  (k)  (  Y(k+1)  -  0*(k)Y(k)  J  Y(k)TP(k)  (16) 
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and 


I  -  j8  (k)Y(k)Y(k)TP(k)        . 


2 
The  above  relationships  identify  an  N  by  N  0  *  matrix  of  N     elements. 

For  the  scalar  observable  case  N(N-l)  of  these  elements  are  already 

2 
known  to  be  zeros  or  ones  .    Hence  N     elements  are  identified  to  learn 

N  unknowns  in  the  N      row.    Thus  a  more  efficient  method  of  the  fore- 
going has  been  derived  by  Yu-Chi  Ho  [  5  ]    and  R.C.K.  Lee  [  4,5  ] 
which  recursively  estimates  the  elements  of  the  N      row. 

0  (k+1)  =  0A(k)  +  P(k)Y(k)  (  z(k+l)  -  0  (k)TY(k)J  j8  (k)  (18) 

P(k)    =  P(k+1)  (l  -  P  (k-l)Y(k-l)Y(k-l)TP(k-l)J  (19) 

where  z(k+l)  =  the  scalar  measurement  at  stage  (k+1) 

a  A       T  th 

0  (k)       =  (Nxl)  column  such  that  0  (k)     is  equal  to  the  N 

row  of  0  *(k) 

Initialization  of  the  recursive  equations  might  be  accomplished  for 
a  -1 

0*(O)  =  S_    S  or  by  using  a  small  sample  set  of  2N+1  or  more 

—  1  —  1 

measurements  and  forming  R(l)  and  R(0)       where  by  0*(O)  =  R(1)R(0) 
and  P(0)  =  R(0)      .     Both  of  these  methods  require  the  complicated  matrix 
inversion  routine.    From  Lee's  work  [4  ]  the  convergence  rate  appears 

A  A 

to  be  satisfactory  for  any  reasonable  estimate  of  0*(O)  or  0  (0).    Ex- 
perience seems  to  bear  this  out  for  the  simulation  model.     The  P  matrix 

A 

ought  to  be  inversely  proportional  to  time  approaching  zero  as  0 
approaches  0  .     P  has  its  most  significant  effect  on  the  initial  rate  of 
convergence.    Experience  indicates  that  if  P  is  initialized  with  a 

c 

reasonably  large  set  of  values  on  the  main  diagonal  (~10  ),  con- 
vergence is  satisfactory.    Additional  work  in  this  area  would  be  useful 
to  define  an  optimal  P(0),  perhaps  in  terms  of  initial  convergence  rate, 
measurement  noise  statistics,  and  identification  error  as  a  function  of 
time. 

Convergence  tests  were  run  using  the  scalar  observable  roll  angles 
generated  by  the  fourth  order  model.    A  ten  second  period,  damping 
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coefficient  of  0.7  and  an  excitation  constant  of  twenty  were  used.    The 
sample  period  was  one  second.    Three  cases  were  investigated  using 
the  same  sample  set. 

The  first  case  was  with  T  *  suppressed  to  zeros  in  all  but  the  last 
element  to  eliminate  numerator  dynamics.    The  convergence  of  600 
samples  is  s.hown  in  fig.  3a,b,c,d  for  each  of  the  fourth  row  elements 
of  0  *.    The  second  case  was  the  same  as  the  first  but  with  Gaussian 
additive  measurement  noise  having  a  standard  deviation  of  one.    Ex- 
cellent convergence  is  shown  in  fig.  4a,b,c,d.     The  third  case  used 
the  entire  T  *  vector  and  no  measurement  noise.    The  residues  are 
apparent  in  fig.  5a,b,c,d.    Comparing  the  eigenvalues  of  0  *  with 

A 

those  of  0  *  with  residues  indicates  a  considerable  shift,  fig.  6.  The 
responses  to  different  sample  rates  in  the  following  section  indicates 
that  the  residues  may  be  minimized  sufficiently  to  provide  acceptable 
root  locations . 

The  true  computed  values  for  the  plant  with  a  10  sec  period  sam- 
pled at  1  sec  intervals  are  as  follows, 

0  CI  0  0 

0  0  10 

0  0  0  1 

-.1722        .9633       -2.177       2.322 


0 


*  = 


r  *  = 


" .105' 

.245 

.500 

.877 

r  *  suppressed  = 


0 

0 

0 

.877 
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Fig.  3.     Recursive  Convergence  of  0  *  vs  Sample  Number,  No  Measurement 
Noise,   10  sec  Dominant  Period  Plant  with  Zero  Initial  Conditions. 
Discrete  Gaussian  Excitation,  <y=  20,  Applied  at  the  Sample 
Interval  of  1  sec .  T  *  Suppressed  (T  *  1,2,3  =  0),  r  =  0 . 7. 
The  Straight  Line  Represents  the  True  Value  of  the  Particular 
0  *  Element,    (a)  $  * (4 , 1)  vs  Sample  Number .   (b)^*(4,2)vs 
Sample  Number,    (c)  0  *(4,3)  vs  Sample  Number,    (d)  0  *(4,4) 
vs  Sample  Number. 

(a)  $  *(4,1)  vs  Sample  Number,    (0  *(4,  1)  =  -0.17217) 
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(b)  #*(4,2)  vs  Sample  Number,   (0  *(4,2)  =  0.96328) 
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(c)  0*(4,3)  vs  Sample  Number,   (0  *(4,3)  =  -2.1772) 
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(d)  0*(4,4)  vs  Sample  Number,   (0  *(4,4)  =  2.3215) 
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Fig.  4.     Recursive  Convergence  of  0  *  vs  Sample  Number  with  Additive 

Gaussian  Measurement  Noise,  cr  =  1 ,  using  the  10  sec  Dominant 
Period  Plant,  Zero  Initial  Conditions.     Sampled  and  Excited  at 
1  sec  Intervals,  T  *  Suppressed  (  T*l  ,2 , 3  =  0),  f  =  0.7.    The 
Straight  Line  Represents  the  True  Value  of  the  Particular  g6  * 
Element,     (a)  %  *(4,1)  vs  Sample  Number,   (b)  <?  *(4,2)  vs  Sample 
Number,   (c)  0  *(4,3)  vs  Sample  Number,    (d)  <fl*(4,4)  vs  Sample 


Number. 


(a)  3  *(4,1)  vs  Sample  Number,   (0*(4,1)  =  -0.17217) 
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(b)  0  *(4,2)  vs  Sample  Number,    (0  *(4,2)  =  0.96328) 
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(c)  £*(4,3)  vs  Sample  Number,    (0  *(4,3)  =  -2.1772) 
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(d)  0  *(4,4)  vs  Sample  Number,   (0  *(4,4)  =  2.3215) 
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Fig.  5.    Recursive  Convergence  of  0  *  vs  Sample  Number,  No  Noise 
using  the  10  sec  Dominant  Period  Plant,  Zero  Initial  Con- 
ditions.    Sampled  and  Excited  at  1  sec  Intervals,  f  =  0.7. 
The  Straight  Line  Represents  the  True  Value  of  the  Particular 
0  *  Element.     The  Full  T  *  vector  was  used  to  Demonstrate 


the  Residues  from  Numerator  Dynamics. 
Number,   (b)  0  *(4,2)  vs  Sample  Number. 
Number,    (d)  0  *(4,4)  vs  Sample  Number. 


(a)  0  *(4,1)  vs  Sample 
(c)  t  *(4,3)  vs  Sample 


(a)  0  *(4,1)  vs  Sample  Number,   (0  *(4,1)  =  -0.17217) 
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(b)  0  *(4,2)  vs  Sample  Number,   (0  *(4,2)  =  0.96328) 
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(c)  A  *(4,3)  vs  Sample  Number,    (0  *(4,3)  =  -2.1772) 
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(d)  0  *(4,4)  vs  Sample  Number,    (0  *(4,4)  =  2.3215) 
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Fig.   6 


Eigenvalues  for  0  * ,  %  *  with  1*1,2,3  =  0,  and  <£  * 
using  the  Complete  T  *  Resulting  in  Residues  for  the 
10  sec  Dominant  Period  Plant  Sampled  and  Excited  at 
1  sec  Intervals.     No  Measurement  Noise. 
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SAMPLE  RATE 
The  choice  of  sample  rate  for  discrete  stochastic  simulation  does 
not  necessarily  follow  the  supposition  that  the  more  rapidly  a  system 
is  sampled  the  "better".    In  some  cases  the  existence  of  an  optimal 
sample  rate  has  been  shown  [  6  ]  .    To  investigate  the  accuracy  of 
identification  by  batch  processing,  a  sample  set  of  500  was  generated 
for  sample  intervals  from  0.02  to  4.0  seconds  in  increments  of  0.02 
seconds.     The  model  was  excited  by  a  different  excitation  set  for  each 

sample  rate  but  with  the  same  statistical  properties,  normal  (0,1) 

th  A 

"white"  noise.     The  N      row  elements  of  0  *  for  each  sample  set  were 

subtracted  from  the  true  0  *  values  and  plotted  versus  sample  inter- 
val, fig.   7a,b,c,d,     For  sample  intervals  less  than  0.06  seconds 
R(0)  went  singular.    Each  element  appears  to  have  its  own  best  sam- 
ple rate  ranging  from  1  to  1.5  seconds  for  a  dominant  10  second  plant. 
Hence  a  rule  of  thumb:    Sample  at  10  times  the  dominant  frequency! 
To  eliminate  the  effect  of  varying  the  excitation  set,  the  above 
procedure  was  repeated  using  the  same  excitation  set  for  each  sample 
period,  0.05  to  10  seconds  in  increments  of  0.05  seconds.    The  re- 

A 

suiting  error  curves  were  the  same  as  illustrated  by  the  0*(4,2)  curve, 
fig.  8a.    The  apparent  stabilization  of  the  error  as  the  sample  period 
approaches  the  natural  period  of  the  system  is  misleading.    The  ab- 
solute values  of  0  *  become  small  due  to  the  large  damping  coefficient 
but  percent  errors  actually  increase. 

It  is  interesting  to  note  that  as  the  dominant  period  is  changed 
the  shape  of  the  error  curve  is  the  same  but  shifted  in  the  direction 
of  change  of  the  period.     The  ratios  of  sample  period  to  natural  period 
remain  constant  at  the  zero  error  intercept  points  .    The  error  curves 

A 

for  the  0*(4,2)  element  for  natural  periods  of  8  and  6  seconds  are 
shown  in  fig.   8b, c.    All  of  these  curves  were  generated  using  the  full 
F  *  vector  and  hence  contain  numerator  dynamics  and  their  inherent 
residues.    Thus  it  is  pertinent  that  there  is  some  optimal  sample  rate 
that  minimizes  residue  errors. 
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Errors  are  not  due  solely  to  numerator  dynamics .    The  original 
contention  that  the  choice  of  sample  period  is  significant  may  be  em- 
phasized by  repeating  the  above  runs  for  the  ten  second  plant  with 
suppressed  T  * .    The  resulting  identification  errors  are  plotted  versus 
sample  period  for  each  element,  fig.   9a,b,c,d.     Errors  may  be  sig- 
nificantly different  if  the  excitation,  W    ,  drives  the  system  at  a 
different  rate  than  the  sampling  (observation)  frequency. 
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Fig.  8.    Error  =  <A  *(4,2)  -  $  *(4,2)  vs  Sample  Period  (0.05  to  10.0  sec) 
using  Batch  Processing  of  500  Samples  from  the  10,  8,  and  6 
sec  Dominant  Period  Plants.     One  Sequence  of  500  Discrete 
Gaussian  Values,  a  =  20,  was  used  as  Excitation  for  all  Sample 
Periods  and  all  3  Plants.    T  *  was  not  Suppressed,   (a)  Error  of 
$  *(4,2)  vs  Sample  Period,   10  sec  Plant,   (b)  Error  of  $* (4,  2)  vs 
Sample  Period,  8  sec  Plant,    (c)  Error  of  $*(A,2)  vs  Sample 
Period,   6  sec  Plant. 

(a)  Error  of  $  *(4,2)  vs  Sample  Period  (sec),    10  sec  Plant. 


TQ 


(b)  Error  of  <t>  *(4,2)  vs  Sample  Period  (sec),   8  sec  Plant. 
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(c)  Error  of  %  *(4,2)  vs  Sample  Period  (sec),   6  sec  Plant. 
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Fig.  9.    Error  =  $  *  -  0  *  vs  Sample  Period  (0.05  to  10.0  sec)  using 
Batch  Processing  of  5  00  Samples  from  the  10  sec  Dominant 
Period  Plant  with  T  *  Suppressed  (T  *  1,2,3  =  0).    Discrete 
Gaussian  Excitation,  a  =  20,  was  applied  at  the  Sample  Rate. 
The  Same  Excitation  Set  was  used  for  each  Different  Sample 
Period,    (a)  Error  of  0  *(4,1)  vs  Sample  Period,   (b)  Error  of 
0  *(4,2)  vs  Sample  Period,   (c)  Error  of  <?*(4,3)  vs  Sample  Period, 
(d)  Error  of  0  *(4,4)  vs  Sample  Period. 

(a)  Error  of  $  *  (4  , 1 )  vs  Sample  Period ,  T  *  i ;  2  ,  3  =  0 
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(b)  Error  of  0  *(4,2)  vs  Sample  Period,  r*if2,3  =  0 
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(c)  Error  of  0  *(4,3)  vs  Sample  Period,  r  *i ,2,3  =  0 
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(d)  Error  of  0  * (4, 4)  vs  Sample  Period,  r  *i,2,3  =  °- 
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APPLICATION  TO  SHIP  MOTION  ESTIMATION 

As  a  result  of  the  foregoing  work  and  associated  simulation,  the 
following  proposal  is  offered  as  a  practical  procedure  for  ship  motion 
estimation. 

Initial  identification  of  the  order  of  the  pitch  and  roll  plants  and 
the  corresponding  0  *'s  is  to  be  evaluated  for  the  free  dynamic  system 
with  the  ship  in  calm  water.    Initial  angle  may  be  applied  with  weights 
or  other  means.    The  data  thus  obtained  may  be  used  as  a  reference  and 
to  initialize  the  stochastic  routines.    The  yaw  plant  may  not  be  iden- 
tified as  a  free  dynamic  system  since  there  is  no  suitable  way  to  apply 
a  controlled  initial  condition  which  will  project  on  all  the  eigenvectors 
in  the  absence  of  a  righting  moment. 

Stochastic  identification  is  to  be  applied  whenever  the  ship  is 
underway  at  sea.     Either  batch  processing  or  the  recursion  equations 
may  be  used.    Batch  processing  should  be  used  at  least  periodically 
to  verify  the  system  order  from  R(0)       and  to  identify  Q  if  required. 

The  recursion  equations  are  more  suitable  for  time  sharing  and  do 
not  require  storage  to  accumulate  large  sample  sets  of  data,     Since 
periodic  identification  requirements  are  anticipated  to  follow  system 
dynamics  ,  recursive  identification  may  be  initiated  using  the  previous 

A 

0*  as  the  best  estimate  and  flagging  in  P(0). 

The  actual  estimation  of  ship  attitude  follows  the  identification 
of  the  0  *'s  (and  Q's).    Equation  (9)  may  be  applied  directly  to  predict. 
Discrete  linear  (Kalman)  filtering  is  available. 

The  development  throughout  this  paper  has  assumed  that  the  only 
reliable  or  convenient  measurements  available  for  identification  were 
discrete  pitch,  roll,  and  yaw  angles.    If  the  number  of  observables 
is  increased,  the  accuracy  of  the  identifications  and  the  estimations 
may  be  expected  to  improve  [  7  ]  . 
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CONCLUSIONS 
Identification  of  ship  motion  dynamics  is  feasible  by  stochastic 
methods . 

The  accuracy  of  the  identification  is  dependent  upon  the  numerator 
dynamics  of  the  true  system  z-transfer  function.    The  contention 
of  Lee  [  4  ]   that  decorrelation  would  improve  accuracy  was  not 
true  for  this  work  or  that  of  Blackner  [  7  ]  . 

The  accuracy  is  dependent  upon  the  sample  period.    An  optimal 
sample  rate  exists. 

Identification  of  an  unexcited  plant  (the  free  dynamic  system) 
having  a  single  initial  condition  is  independent  of  r*  and  hence 
is  not  affected  by  numerator  dynamics.    This  should  be  the  source 
of  the  most  reliable  identification  as  a  reference. 
Further  work  should  be  done  using  actual  ship  motion  data  to 
complete  the  identification  evaluation  and  to  investigate  the 
performance  of  discrete  linear  (Kalman)  filtering  as  the  final 
stage  of  ship  motion  estimation. 
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APPENDIX 
CDC  FORTRAN  63  Programs  and  Subroutines  for  Digital  Sim- 
ulation and  Identification. 

PHITEST:  Main  program  generating  the  simulation  model 

and  controlling  the  identification. 
PHIDENT:         Recursive  identification  routine  for  the  scalar 

observable  case. 
PHICALL:  Batch  processing  for  general  observable  case. 

Identifies  0  *,  N,  and  Q.   (GAUSS3  is  identical 

to  RECIP) 
PHIDEL:  Computes  standard  form  0  and  T  matrices  for 

sample  data  systems. 
PHIDEL1:  Computes  0  *  and  T  *  matrices  for  sample 

data  systems . 
PSD:  Computes  Power  Spectral  Density  from  auto- 

correlation functions . 
RECIP:  Matrix  inversion  routine. 

ADD,  SUB,  TRANS,  PROD,   PRINTER,  and  MATREAD: 

Convenient  subroutines  for  the  indicated 

matrix  operations. 
PTPLOT:  Graphs  and  prints  out  2  dimensional  data  on 

line  printer. 
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